library(tidyverse)
library(modelsummary)
library(MASS)
library(jtools)


d <- read.csv("broome_data.csv")

# vectors of variables needed for commands below
covars <- paste(c("female", "pid7",  "age", 
                  "educ", "race_nonwhite", "interest", 
                  "turnout_likelihood", "voteintention", 
                  "percentdemvoter", "bigmoney_importance", "obama_approve", 
                  "fracking_support", "limitcontributions"), collapse = " + ")
cv <- c("female", "pid7",  "age", 
        "educ", "race_nonwhite", "interest", 
        "turnout_likelihood", "voteintention", 
        "percentdemvoter", "bigmoney_importance", "obama_approve", 
        "fracking_support", "limitcontributions")


##################
#### analysis ####
##################



#### Table 2 ####
## Amendment policy attitudes ##
## wave 2 #

round(prop.table(table(d$treatment, d$cf_amend_w2), margin = 1), 4)*100

# Ordered Logit
m1 <- polr(as.formula(factor(cf_amend_w2) ~ treatment), 
           data = d, Hess=TRUE)

m2 <- polr(as.formula(paste0("factor(cf_amend_w2) ~ treatment +", covars)), 
           data = d, Hess=TRUE)

## wave 3 ##

round(prop.table(table(d$treatment, d$cf_amend_w3), margin = 1), 4)*100

# Ordered Logit
m3 <- polr(as.formula(factor(cf_amend_w3) ~ treatment), 
           data = d, Hess=TRUE)

m4 <- polr(as.formula(paste0("factor(cf_amend_w3) ~ treatment +", covars)), 
           data = d, Hess=TRUE)

modelsummary(list(m1, m2, m3, m4), vcov = ~cluster_id, stars = T,  
             coef_map = c("treatmentAnti-Amendment" = "Anti-Amendment",
                          "treatmentPro-Amendment" = "Pro-Amendment",
                          "female" = "Female",
                          "pid7" = "Party ID",
                          "age" = "Age",
                          "educ" = "Education",
                          "race_nonwhite" = "Nonwhite",
                          "interest" = "Political Interest",
                          "turnout_likelihood" = "Turnout Likelihood",
                          "percentdemvoter" = "% Dem. Voter",
                          "obama_approve" = "Obama Approval",
                          "fracking_support" = "Fracking Support",
                          "limitcontributions" = "Limit Contributions",
                          "0|1" = "0|1 cutoff",
                          "1|2" = "1|2 cutoff",
                          "2|3" = "2|3 cutoff",
                          "3|4" = "3|4 cutoff"),
             notes = "Ordered logit regression estimates and robust standard errors clustered by household.")

#### Figure 1 ######
## wave 2 ##
d |> mutate(cf_amend_w2_char = case_match(cf_amend_w2,
                                          0 ~ "Strongly oppose",
                                          1 ~ "Somewhat oppose",
                                          2 ~ "Neither oppose \nnor support",
                                          3 ~ "Somewhat support",
                                          4 ~ "Strongly support"),
            cf_amend_w2_char = fct_relevel(cf_amend_w2_char, 
                                           c("Strongly oppose", "Somewhat oppose",
                                             "Neither oppose \nnor support", 
                                             "Somewhat support", 
                                             "Strongly support"))) |>  
      filter(!is.na(cf_amend_w2_char)) |>
      group_by(treatment) |>
      count(cf_amend_w2_char) |>
      reframe(cf_amend_w2_char = cf_amend_w2_char,
              Percent = n/sum(n), 
              se = sqrt(Percent*(1-Percent)/(sum(n)))) |>
      ggplot(aes(cf_amend_w2_char, Percent, fill = treatment)) + 
      geom_col(position = position_dodge()) +
      geom_errorbar(aes(ymin = Percent - 1.96*se,
                        ymax = Percent + 1.96*se), 
                    position = position_dodge(width = .9), 
                    width = .2,
                    alpha = .7) + 
      scale_y_continuous(labels = scales::percent) +
      scale_fill_grey(labels = c("Fracking", "Anti-Reform", "Pro-Reform"), 
                      name = "Treatment") + 
      labs(x = "Support for \'Reasonable\' Limits") + 
      theme_bw()

#ggsave("Figure1.png", dpi = 900, width = 7)

#### Figure 2 ####
## wave 3 ##
d |> mutate(cf_amend_w3_char = case_match(cf_amend_w3,
                                          0 ~ "Strongly oppose",
                                          1 ~ "Somewhat oppose",
                                          2 ~ "Neither oppose \nnor support",
                                          3 ~ "Somewhat support",
                                          4 ~ "Strongly support"),
            cf_amend_w3_char = fct_relevel(cf_amend_w3_char, 
                                    c("Strongly oppose", "Somewhat oppose",
                                    "Neither oppose \nnor support", 
                                 "Somewhat support", 
                                 "Strongly support"))) |>  
      filter(!is.na(cf_amend_w3_char)) |>
      group_by(treatment) |>
      count(cf_amend_w3_char) |>
      reframe(cf_amend_w3_char = cf_amend_w3_char,
              Percent = n/sum(n), 
              se = sqrt(Percent*(1-Percent)/(sum(n)))) |>
      ggplot(aes(cf_amend_w3_char, Percent, fill = treatment)) + 
      geom_col(position = position_dodge()) +
      geom_errorbar(aes(ymin = Percent - 1.96*se,
                        ymax = Percent + 1.96*se), 
                    position = position_dodge(width = .9), 
                    width = .2,
                    alpha = .7) + 
      scale_y_continuous(labels = scales::percent) +
      scale_fill_grey(labels = c("Fracking", "Anti-Reform", "Pro-Reform"), 
                      name = "Treatment") + 
      labs(x = "Support for Overturning Citizens United") + 
      theme_bw()

#ggsave("Figure2.png", dpi = 900, width = 7)      




#### Figure 3 ####
## wave 2 ##

m2 <- polr(as.formula(paste0("factor(cf_imp2_w2) ~ treatment +", 
                             covars)), data = d, Hess = T) 

# placebo
newd <- mutate(m2$model, 
               female = 1, 
               treatment = "0Fracking",
               pid7 = mean(m2$model$pid7), 
               age = mean(m2$model$age),
               educ = mean(m2$model$educ),
               race_nonwhite = 0,
               interest = mean(m2$model$interest),
               turnout_likelihood = mean(m2$model$turnout_likelihood),
               voteintention = mean(m2$model$voteintention), 
               percentdemvoter = mean(m2$model$percentdemvoter),
               bigmoney_importance = mean(m2$model$bigmoney_importance),
               obama_approve = mean(m2$model$obama_approve),
               fracking_support = mean(m2$model$fracking_support),
               limitcontributions = mean(m2$model$limitcontributions))

yhat <- predict(m2, newdata = newd, type = "probs")
fitted0 <- data.frame(fit = colMeans(yhat))
fitted0$condition = "Placebo"
fitted0$label <- c("Not at all", "Not very", "Somewhat", "Very")


# pro-amendment
newd <- mutate(m2$model, 
               female = 1, 
               treatment = "Pro-Amendment",
               pid7 = mean(m2$model$pid7), 
               age = mean(m2$model$age),
               educ = mean(m2$model$educ),
               race_nonwhite = 0,
               interest = mean(m2$model$interest),
               turnout_likelihood = mean(m2$model$turnout_likelihood),
               voteintention = mean(m2$model$voteintention), 
               percentdemvoter = mean(m2$model$percentdemvoter),
               bigmoney_importance = mean(m2$model$bigmoney_importance),
               obama_approve = mean(m2$model$obama_approve),
               fracking_support = mean(m2$model$fracking_support),
               limitcontributions = mean(m2$model$limitcontributions))

yhat <- predict(m2, newdata = newd, type = "probs")
colMeans(yhat)
fitted1 <- data.frame(fit = colMeans(yhat))
fitted1$condition = "Pro-Amendment"
fitted1$label <- c("Not at all", "Not very", "Somewhat", "Very")

# anti-amendment
newd <- mutate(m2$model, 
               female = 1, 
               treatment = "Anti-Amendment",
               pid7 = mean(m2$model$pid7), 
               age = mean(m2$model$age),
               educ = mean(m2$model$educ),
               race_nonwhite = 0,
               interest = mean(m2$model$interest),
               turnout_likelihood = mean(m2$model$turnout_likelihood),
               voteintention = mean(m2$model$voteintention), 
               percentdemvoter = mean(m2$model$percentdemvoter),
               bigmoney_importance = mean(m2$model$bigmoney_importance),
               obama_approve = mean(m2$model$obama_approve),
               fracking_support = mean(m2$model$fracking_support),
               limitcontributions = mean(m2$model$limitcontributions))

yhat <- predict(m2, newdata = newd, type = "probs")
colMeans(yhat)
fitted2 <- data.frame(fit = colMeans(yhat))
fitted2$condition = "Anti-Amendment"
fitted2$label <- c("Not at all", "Not very", "Somewhat", "Very")

fitted <- bind_rows(fitted0, fitted1, fitted2)

fitted |> mutate(condition = fct_relevel(condition, 
                                         c("Placebo", "Anti-Amendment")),
                 se = sqrt(fit*(1-fit)/(nobs(m2)))) |>
      ggplot(aes(label, fit, fill = condition)) + 
      geom_col(position = position_dodge()) +
      geom_errorbar(aes(ymin = fit - 1.96*se,
                       ymax = fit + 1.96*se), 
                    position = position_dodge(width = .9), 
                    width = .2,
                    alpha = .7) + 
      scale_y_continuous(labels = scales::percent) +
      scale_fill_grey(labels = c("Fracking", "Anti-Reform", "Pro-Reform"), 
                      name = "Treatment") + 
      labs(x = "Personal Importance of Campaign Finance",
           y = "Predicted Probability (in %)") + 
      theme_bw()


#### Table 3 ####
## wave 3 ##
# Ordered Logit
m1 <- polr(as.formula(factor(cf_imp2_w2) ~ treatment), data = d, Hess = T)

m2 <- polr(as.formula(paste0("factor(cf_imp2_w2) ~ treatment +", 
                             covars)), data = d, Hess = T) 

m3 <- polr(as.formula(factor(cf_imp2_w3) ~ treatment), data = d, Hess = T)

m4 <- polr(as.formula(paste0("factor(cf_imp2_w3) ~ treatment +", 
                             covars)), data = d, Hess = T) 

modelsummary(list(m1, m2, m3, m4), vcov = ~cluster_id, stars = T, 
             coef_map = c("treatmentAnti-Amendment" = "Anti-Amendment",
                          "treatmentPro-Amendment" = "Pro-Amendment",
                          "female" = "Female",
                          "pid7" = "Party ID",
                          "age" = "Age",
                          "educ" = "Education",
                          "race_nonwhite" = "Nonwhite",
                          "interest" = "Political Interest",
                          "turnout_likelihood" = "Turnout Likelihood",
                          "percentdemvoter" = "% Dem. Voter",
                          "obama_approve" = "Obama Approval",
                          "fracking_support" = "Fracking Support",
                          "limitcontributions" = "Limit Contributions",
                          "0|1" = "0|1 cutoff",
                          "1|2" = "1|2 cutoff",
                          "2|3" = "2|3 cutoff",
                          "3|4" = "3|4 cutoff"))

#### Table 4 #### 
## prior-attitude moderation ##
# anti-and pro-limit proportions
prop.table(table(d$limitcontributions))

## wave 2 ##
m1 <- polr(as.formula(paste0("factor(cf_imp2_w2) ~ treatment * limitcontributions")), 
           data = d, Hess = T)

m2 <- polr(as.formula(paste0("factor(cf_imp2_w2) ~ treatment * limitcontributions +", 
                             covars)), 
           data = d, Hess = T)

## wave 3 ##
m3 <- polr(as.formula(paste0("factor(cf_imp2_w3) ~ treatment * limitcontributions")), 
           data = d, Hess = T)

m4 <- polr(as.formula(paste0("factor(cf_imp2_w3) ~ treatment * limitcontributions +", 
                             covars)), 
           data = d, Hess = T)


modelsummary(list(m2, m4), vcov = ~cluster_id, stars = T, 
             coef_map = c("treatmentAnti-Amendment" = "Anti-Amendment x Anti-Limit",
                          "treatmentAnti-Amendment:limitcontributions" = "Anti-Amendment x Pro-Limit",
                          "treatmentPro-Amendment" = "Pro-Amendment x Anti-Limit",
                          "treatmentPro-Amendment:limitcontributions" = "Pro-Amendment x Pro-Limit",
                          "female" = "Female",
                          "pid7" = "Party ID",
                          "age" = "Age",
                          "educ" = "Education",
                          "race_nonwhite" = "Nonwhite",
                          "interest" = "Political Interest",
                          "turnout_likelihood" = "Turnout Likelihood",
                          "percentdemvoter" = "% Dem. Voter",
                          "obama_approve" = "Obama Approval",
                          "fracking_support" = "Fracking Support",
                          "limitcontributions" = "Limit Contributions",
                          "0|1" = "0|1 cutoff",
                          "1|2" = "1|2 cutoff",
                          "2|3" = "2|3 cutoff",
                          "3|4" = "3|4 cutoff"))




#### Table 5 ####
m4 <- polr(as.formula(paste0("factor(cf_amend_w2) ~ treatment * limitcontributions +", covars)), 
           data = d, Hess=TRUE)

m6 <- polr(as.formula(paste0("factor(cf_amend_w3) ~ treatment * limitcontributions +", covars)), 
           data = d, Hess=TRUE)

modelsummary(list(m4, m6), vcov = ~cluster_id, stars = T, 
             coef_map = c("treatmentAnti-Amendment" = "Anti-Amendment x Anti-Limit",
                          "treatmentAnti-Amendment:limitcontributions" = "Anti-Amendment x Pro-Limit",
                          "treatmentPro-Amendment" = "Pro-Amendment x Anti-Limit",
                          "treatmentPro-Amendment:limitcontributions" = "Pro-Amendment x Pro-Limit",
                          "female" = "Female",
                          "pid7" = "Party ID",
                          "age" = "Age",
                          "educ" = "Education",
                          "race_nonwhite" = "Nonwhite",
                          "interest" = "Political Interest",
                          "turnout_likelihood" = "Turnout Likelihood",
                          "percentdemvoter" = "% Dem. Voter",
                          "obama_approve" = "Obama Approval",
                          "fracking_support" = "Fracking Support",
                          "limitcontributions" = "Limit Contributions",
                          "0|1" = "0|1 cutoff",
                          "1|2" = "1|2 cutoff",
                          "2|3" = "2|3 cutoff",
                          "3|4" = "3|4 cutoff"))


#### Figure 4 ####
## wave 3 ##
m4 <- polr(as.formula(paste0("factor(cf_imp2_w3) ~ treatment * limitcontributions +", 
                             covars)), 
           data = d, Hess = T)

# placebo; limit = 0
newd <- mutate(m4$model, 
               female = 1, 
               treatment = "0Fracking",
               pid7 = mean(m4$model$pid7), 
               age = mean(m4$model$age),
               educ = mean(m4$model$educ),
               race_nonwhite = 0,
               interest = mean(m4$model$interest),
               turnout_likelihood = mean(m4$model$turnout_likelihood),
               voteintention = mean(m4$model$voteintention), 
               percentdemvoter = mean(m4$model$percentdemvoter),
               bigmoney_importance = mean(m4$model$bigmoney_importance),
               obama_approve = mean(m4$model$obama_approve),
               fracking_support = mean(m4$model$fracking_support),
               limitcontributions = 0)

yhat <- predict(m4, newdata = newd, type = "probs")
fitted0 <- data.frame(fit = colMeans(yhat))
fitted0$condition = "Placebo"
fitted0$prior = "Anti-Limit"
fitted0$label <- c("Not at all", "Not very", "Somewhat", "Very")

# placebo; limit = 1
newd <- mutate(m4$model, 
               female = 1, 
               treatment = "0Fracking",
               pid7 = mean(m4$model$pid7), 
               age = mean(m4$model$age),
               educ = mean(m4$model$educ),
               race_nonwhite = 0,
               interest = mean(m4$model$interest),
               turnout_likelihood = mean(m4$model$turnout_likelihood),
               voteintention = mean(m4$model$voteintention), 
               percentdemvoter = mean(m4$model$percentdemvoter),
               bigmoney_importance = mean(m4$model$bigmoney_importance),
               obama_approve = mean(m4$model$obama_approve),
               fracking_support = mean(m4$model$fracking_support),
               limitcontributions = 1)

yhat <- predict(m4, newdata = newd, type = "probs")
fitted01 <- data.frame(fit = colMeans(yhat))
fitted01$condition = "Placebo"
fitted01$prior = "Pro-Limit"
fitted01$label <- c("Not at all", "Not very", "Somewhat", "Very")


# pro-amendment; limit = 0
newd <- mutate(m4$model, 
               female = 1, 
               treatment = "Pro-Amendment",
               pid7 = mean(m4$model$pid7), 
               age = mean(m4$model$age),
               educ = mean(m4$model$educ),
               race_nonwhite = 0,
               interest = mean(m4$model$interest),
               turnout_likelihood = mean(m4$model$turnout_likelihood),
               voteintention = mean(m4$model$voteintention), 
               percentdemvoter = mean(m4$model$percentdemvoter),
               bigmoney_importance = mean(m4$model$bigmoney_importance),
               obama_approve = mean(m4$model$obama_approve),
               fracking_support = mean(m4$model$fracking_support),
               limitcontributions = 0)

yhat <- predict(m4, newdata = newd, type = "probs")
colMeans(yhat)
fitted1 <- data.frame(fit = colMeans(yhat))
fitted1$condition = "Pro-Amendment"
fitted1$prior = "Anti-Limit"
fitted1$label <- c("Not at all", "Not very", "Somewhat", "Very")

# pro-amendment; limit = 1
newd <- mutate(m4$model, 
               female = 1, 
               treatment = "Pro-Amendment",
               pid7 = mean(m4$model$pid7), 
               age = mean(m4$model$age),
               educ = mean(m4$model$educ),
               race_nonwhite = 0,
               interest = mean(m4$model$interest),
               turnout_likelihood = mean(m4$model$turnout_likelihood),
               voteintention = mean(m4$model$voteintention), 
               percentdemvoter = mean(m4$model$percentdemvoter),
               bigmoney_importance = mean(m4$model$bigmoney_importance),
               obama_approve = mean(m4$model$obama_approve),
               fracking_support = mean(m4$model$fracking_support),
               limitcontributions = 1)

yhat <- predict(m4, newdata = newd, type = "probs")
colMeans(yhat)
fitted2 <- data.frame(fit = colMeans(yhat))
fitted2$condition = "Pro-Amendment"
fitted2$prior = "Pro-Limit"
fitted2$label <- c("Not at all", "Not very", "Somewhat", "Very")


# anti-amendment; limit = 0
newd <- mutate(m4$model, 
               female = 1, 
               treatment = "Anti-Amendment",
               pid7 = mean(m4$model$pid7), 
               age = mean(m4$model$age),
               educ = mean(m4$model$educ),
               race_nonwhite = 0,
               interest = mean(m4$model$interest),
               turnout_likelihood = mean(m4$model$turnout_likelihood),
               voteintention = mean(m4$model$voteintention), 
               percentdemvoter = mean(m4$model$percentdemvoter),
               bigmoney_importance = mean(m4$model$bigmoney_importance),
               obama_approve = mean(m4$model$obama_approve),
               fracking_support = mean(m4$model$fracking_support),
               limitcontributions = 0)

yhat <- predict(m4, newdata = newd, type = "probs")
colMeans(yhat)
fitted3 <- data.frame(fit = colMeans(yhat))
fitted3$condition = "Anti-Amendment"
fitted3$prior = "Anti-Limit"
fitted3$label <- c("Not at all", "Not very", "Somewhat", "Very")

# anti-amendment; limit = 1
newd <- mutate(m4$model, 
               female = 1, 
               treatment = "Anti-Amendment",
               pid7 = mean(m4$model$pid7), 
               age = mean(m4$model$age),
               educ = mean(m4$model$educ),
               race_nonwhite = 0,
               interest = mean(m4$model$interest),
               turnout_likelihood = mean(m4$model$turnout_likelihood),
               voteintention = mean(m4$model$voteintention), 
               percentdemvoter = mean(m4$model$percentdemvoter),
               bigmoney_importance = mean(m4$model$bigmoney_importance),
               obama_approve = mean(m4$model$obama_approve),
               fracking_support = mean(m4$model$fracking_support),
               limitcontributions = 1)

yhat <- predict(m4, newdata = newd, type = "probs")
colMeans(yhat)
fitted4 <- data.frame(fit = colMeans(yhat))
fitted4$condition = "Anti-Amendment"
fitted4$prior = "Pro-Limit"
fitted4$label <- c("Not at all", "Not very", "Somewhat", "Very")


fitted <- bind_rows(fitted0, fitted01, fitted1, fitted2, fitted3, fitted4)

fitted |> mutate(condition = fct_relevel(condition, 
                                         c("Placebo", "Anti-Amendment")),
                 se = sqrt(fit*(1-fit)/(nobs(m4)))) |>
      ggplot(aes(label, fit, fill = condition)) + 
      geom_col(position = position_dodge()) +
      geom_errorbar(aes(ymin = fit - 1.96*se,
                        ymax = fit + 1.96*se), 
                    position = position_dodge(width = .9), 
                    width = .2,
                    alpha = .7) + 
      scale_y_continuous(labels = scales::percent) +
      scale_fill_grey(labels = c("Fracking", "Anti-Reform", "Pro-Reform"), 
                      name = "Treatment") + 
      labs(x = "Personal Importance of Campaign Finance",
           y = "Predicted Probability (in %)") + 
      theme_bw() +
      facet_grid(~prior)


#### balance table ####
foo <- d |> mutate(wave = as.factor(case_when(
      wave2 == 1 ~ 2,
      wave3 == 1 ~ 3)),
      bigmoney_1 = as.numeric(bigmoney_importance == 4)) |>
      dplyr::select(all_of(cv), treatment, wave, bigmoney_1) 



datasummary(female + pid7 + age + educ + race_nonwhite + interest + 
                  turnout_likelihood + voteintention + percentdemvoter + 
                  bigmoney_1 + obama_approve + fracking_support + 
                  limitcontributions ~ wave * treatment *(mean), 
                    data = foo)









